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We analyze the slow, glassy structural relaxation as measured through collective and tagged- 
particle density correlation functions obtained from Brownian dynamics simulations for a poly- 
disperse system of quasi-hard spheres in the framework of the mode-coupling theory of the glass 
transition (MCT). Asymptotic analyses show good agreement for the collective dynamics when 
polydispersity effects are taken into account in a multi-component calculation, but qualitative dis- 
agreement at small q when the system is treated as effectively monodisperse. The origin of the 
different small-g behaviour is attributed to the interplay between interdiffusion processes and struc- 
tural relaxation. Numerical solutions of the MCT equations are obtained taking properly binned 
partial static structure factors from the simulations as input. Accounting for a shift in the critical 
density, the collective density correlation functions are well described by the theory at all densities 
investigated in the simulations, with quantitative agreement best around the maxima of the static 
structure factor, and worst around its minima. A parameter- free comparison of the tagged-particle 
dynamics however reveals large quantiative errors for small wave numbers that are connected to the 
well-known decoupling of self-diffusion from structural relaxation and to dynamical heterogeneities. 
While deviations from MCT behaviour are clearly seen in the tagged-particle quantities for densi- 
ties close to and on the liquid side of the MCT glass transition, no such deviations are seen in the 
collective dynamics. 

PACS numbers: 64.70.Pf,82.70.Dd,61.20.Lc 



I. INTRODUCTION 

Understanding the slow dynamical processes occur- 
ring in supercooled glass-forming liquids is still one of 
the challenges in condensed matter physics. The mode- 
coupling theory of the glass transition (MCT), intro- 
duced in 1984 by Bengtzelius, Gotze, and Sjolander and 
Leutheusser [l|, Q , provides a quantitative description of 
the initial slowing down of structural relaxation, when 
approaching the glassy state from the liquid region. The 
theory predicts an ideal glass transition whose signature 
is a two-step relaxation of dynamical density correlation 
functions. It arises from a divergence of two time scales 
connected with the intermediate relaxation ascribed to 
relaxation of particles inside their neighbor-cages (the /3 
relaxation), and with the final escape of particles from 
their initial positions that restores liquid-like motion (the 
a relaxation process). Two hallmarks of glassy dynam- 
ics, viz. nonexponential, "stretched-exponential" a relax- 
ation and its scaling (also known as "time-temperature 
superposition principle"), are predicted as asymptotic re- 
sults of MCT. The theory stimulated many experiments 
specifically to address the dynamical window for which 
the theory was designed; among them dynamic light scat- 
tering performed on colloidal systems, Brillouin and neu- 
tron scattering, dielectric spectroscopy, and computer 
simulation studies (see Ref. [S-Q for reviews). 

However, MCT is based on the ad-hoc assumption that 
the fluctuating forces (the longitudinal projections of the 
microscopic stresses for a particular wave vector q) are 



governed entirely by the dynamics of density pair fluctu- 
ations. To close the equations, one then further approxi- 
mates a dynamical four-point average through a product 
of density correlation functions. Even though the theory 
has had many successes in describing some key features of 
the slow dynamics, often validated through comparison 
of its asymptotic formula with experimental and simu- 
lation data, the accuracy of the MCT approximation is 
still largely unknown. In molecular glass formers, the fact 
that relaxation times do not diverge at the MCT transi- 
tion, but continue to grow smoothly in a regime where 
motion is thought to be no longer liquid-like but gov- 
erned by activated, so-called "hopping" , processes, is the 
most widely criticized manifestation of the approximate 
nature of MCT. Another commonly quoted feature that 
is not contained in the theory are the non-Gaussian dis- 
tributions of particle displacement discussed in terms of 
dynamical heterogeneity 0, Q • This effect is often linked 
to the appearance of a decoupling of viscous and diffusive 
time scales - the breakdown of the Stokes-Einstein rela- 
tion, although it could be argued that not its breakdown 
at low temperatures, but rather its validity at higher tem- 
peratures in complex glass formers is the surprising fea- 
ture. 

If hopping processes are indeed what is missing in 
MCT, checking the feasibility of creating an atomistic 
model system where such effects are absent, is an obvi- 
ous thing to do. Following the pioneering dynamic light- 
scattering experiments on colloidal hard-sphere-like sus- 
pensions by Pusey, van Megen and coworkers (9l-[l3| , hard 
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spheres with Brownian short-time motion have some- 
times been quoted in this regard. Yet, this exceptional 
nature of the hard-sphere glass transition has been chal- 
lenged based on computer simulation of tagged-particle 
density correlation functions IJ]. In this contribution, 
we wish to analyze the situation further by shifting focus 
from the incoherent quantities to the collective density 
correlation functions that are closer to the framework of 
MCT. We find that, a number of commonly discussed 
shortcomings of the theory appears only in the tagged- 
particle but not in the collective dynamics. 

One strength of MCT is that it allows, in principle, to 
predict detailed information on the slow dynamics when 
given only the particles' interaction potentials, in the 
form of the static structure factor, as input. For real- 
world glass formers, generally mixtures or moderately 
complicated organic molecules, resolving all the partial 
static structure information required is a formidable task. 
Even more so if one is interested not only in the static 
structure, but also the corresponding dynamical relax- 
ation functions. Thus, in many experimental studies, 
MCT results were taken either from asymptotic expan- 
sions (that cannot address the molecular details and 
preasymptotic corrections, which may be strong), or from 
schematic simplifications of the theory's equations (re- 
sulting in a set of fit parameters whose physical meaning 
is rather unclear). Only recently has it become possible 
to perform MCT calculations based on actual experimen- 
tally measured partial structure factors, due to advances 
in neutron scattering techniques on liquid metallic melts 

Thus, testing the "full MCT" , that is, putting to test 
the dynamics as calculated within the theory from the 
static structure factor (without invoking asymptotic or 
schematic limits of the theory's equations) against the 
measured one, is a task for molecular dynamics simu- 
lations, and has been performed on the standar d g lass- 
forming binary Lennard- Jones mixture [l^ . [iGl - llSj . on 
hard-sphere mixtures fl9l . [20| , soft spheres with short- 
ranged attraction [2l|, |22| . and in more complicated sys - 
tems such as network- forming stron g li quids psi . [2^ . 
metallic glasses H^, polymer melts |26l|. or computer 
models of organic glass formers such as ortho-terphenyl 
[13, 113. As it turns out, these systems arc already quite 
demanding for MCT, although the theory fares well in 
a qualitative description of the dynamical phenomena, 
sometimes even quantiatively (most notably, Ref. [23j . 
where also static triplet correlation functions have been 
extracted from simulation and fed into MCT, addressing 
a term in the MCT equations whose existence is often 
silently ignored). 

The most simple model for a classical dense liquid is 
arguably the hard sphere system. In a previous study, 
we addressed a test of MCT for this system partially, 
by comparing MCT and molecular dynamics simulations 
for a polydispcrsc quasi-hard sphere system [29| . There, 
however, computational limitations in acquiring the de- 
sired statistics restricted the discussion to the single- 



particle dynamics (in form of the incoherent density cor- 
relation functions and quantities derived from it, such as 
the mean-squared displacements or diffusion coefficients) . 
It should be stressed that MCT is, in its very essence, a 
theory for the collective slowing down caused by a feed- 
back mechanism for the collective, or coherent density 
correlation functions. Calculating tagged-particle dy- 
namics from this viewpoint involves an additional level of 
(MCT-approximate) equations, and can thus be viewed 
as a more indirect way of testing the theory. In this pa- 
per, we complete the task of Ref. [2^ by detailing a com- 
prehensive, quantitative comparison of MCT with molec- 
ular dynamics computer simulations for the same quasi 
hard-sphere system on the level of the collective density 
correlation functions. Additionally, while in Ref. [l^ an 
approximate liquid-state theory for the static structure 
factor input to MCT was used (the Percus-Yevick ap- 
proximation), we avoid this additional non-MCT level 
of approximations by using directly the simulated static 
structure factors. 

Comparing the dynamics to MCT, it should be rec- 
ognized that the theory focuses on the slow structural 
relaxation, mistreating the short-time dynamics as gov- 
erned by uncorrelated binary collisions (whose inclusing 
into the theory is not straightforward [la, [l3l ) . It is thus 
desirable to minimize the influence of this short-time dy- 
namics, in particular since its details do not change those 
of the long-time relaxation [s^, HH ; it turns out that this 
is achieved by including stochastic noise in the simulated 
equations of motion, leading to a Langevin-dynamics 
simulation. Standard molecular-dynamics integration is 
then most easily implemented using a regular soft-sphere 
potential, V{r) cx r~^^ . For such steep power- law poten- 
tials, it is known that the influence of slight 'softness' on 
the dynamics can be mapped to an effective density in 
a mapping that takes into account the according shift of 
the freezing point [s^. In addition, the stochastic dy- 
namics provides a link to experimental data on colloidal 
hard-sphere-like systems, whose Brownian short-time dy- 
namics we mimick in our simulations. 

Monodispcrse hard- or soft-sphere systems beyond the 
freezing point readily crystallize in simulation. To avoid 
this, we take the spheres' radii to be polydisperse, evenly 
sampled from a narrow distribution just wide enough 
to suppress crystallization on the time scales considered 
in our simulations. It is well known that already small 
polydispersities are very efficient in slowing down nucle- 
ation events dramatically [H, [sj] ; polydispersity is also 
inherent to most colloidal suspensions, making it a nat- 
ural feature to consider. Binary mixtures are another 
common way of circumventing unwanted crystallization, 
but polydisperse systems have the advantage that one 
can meaningfully construct species-averaged total corre- 
lation functions that are not too different from the in- 
dividual partial correlations. In principle, this allows 
to greatly simplify the discussion, by applying the origi- 
nal one-component formulation of MCT. However, poly- 
dispersity may play an interesting role when comparing 
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theory and simulations in particular at small wave num- 
bers q. We will discuss these points in detail below, in- 
cluding three- and five-component moment approxima- 
tions to the polydisperse radius distribution in multi- 
component MCT [H, [3^. Comparing with the one- 
component MCT, this addresses the question of static 
versus dynamic averaging: while the true dynamics of the 
system can be mapped onto an effective one-component 
system only by averaging at the level of the dynamical 
correlation functions (a procedure to which we will refer 
as post-averaging) , it is tempting to perform such averag- 
ing over a narrow polydispersity distribution already on 
the level of the static structure factor (pre-averaging) . 
However, as we will discuss, the nonlinear feedback ef- 
fects of the dynamics pose a limit to the validity of such 
an approach. 

The paper is organized as follows: In Section [IT] the 
features needed in the further analysis of both MCT and 
simulation are reported for reference. Section |llT] demon- 
strates purely asymptotic analyses of the simulation data, 
while Sec. IIVI turns to the full MCT description of the 
simulated dynamical correlation functions. Finally, in 
Sec.|V]we summarize. 

II. SIMULATION AND MCT 

A. Molecular-dynamics simulation 

We perform strongly damped molecular-dynamics 
(Langevin-dynamics) simulations mimicking colloidal 
Brownian dynamics (BD). The core-core respulsion be- 
tween particles i and j is given by 

where di2 is the center-to-center distance, dij = (di + 
dj)/2, with di the diameters of the particles, sampled 
from a uniform distribution centered on the mean di- 
ameter d with half- width S = O.ld. Such a soft-sphere 
system has only one control parameter given by a spe- 
cific combination of number density p and temperature 
T, r = [m-3K-i2] [13. We vary T by keeping the 

temperature fixed and changing the system's density. It 
has been shown (32j that the exponent n = 36 of the 
inverse power-law potential is large enough to effectively 
approximate hard-sphere behavior. 

The equation of motion for particle j is given by the 
Langevin equation, 

i 

which contains the direct forces between particles Fij, 
and stochastic and friction forces, with friction coeffi- 
cient 7, modeling interaction with a solvent. Assuming 
Stokes friction, its value would be connected to the sol- 
vent viscosity 77s, 7 = Sird'qs, where d is the hydrody- 
namic diameter of the particle (approximated as equal 



for all particles, since the spread in the di is small). 
The random forces fulfill the fiuctuation-dissipation the- 
orem, (Ij(t)lj(t')) = 6kBTjS{t - t')5ij. Let us note that 
with the value of 7 chosen in our simulations, the short- 
time dynamics visible in the correlators and in the mean- 
squared displacement is not yet completely overdamped, 
i.e., it is not strictly diffusive, but rather strongly damped 
ballistic. Since it is not our aim to investigate the very 
short-time dynamical features of the simulations, this will 
not be discussed in the following. 

Equilibration runs were performed with undamped 
Newtonian dynamics in all cases, since the damping in- 
troduces a slowing down in the overall time scale. N = 
1000 particles are simulated in a cubic box with standard 
periodic boundary conditions. Lengths are measured in 
units of the mean diameter d, the particle mass m = 1, 
and the unit of time is fixed setting the thermal veloc- 
ity to wth = (/cBr/m)^/^ = l/-\/3. The damping was 
chosen as 7 = 20 in these dimensionless units, and the 
equations of motion were integrated following Heun's al- 
gorithm [13 with a time step of 5t = 0.0005. Due to the 
polydispersity, crystallization did not occur in the runs 
that have been analyzed for the following discussions. 
Crystallization was monitored through the orientational 
order parameter Qq [H,!!^] . It was found that 8 out of 18 
nms for ip = 0.58, and 16 out of 26 runs for = 0.585 in 
fact did crystallize and had to be excluded. Density is re- 
ported as volume fraction, = {-k /Q)d^ [l + (5^] p, where 
the polydispersity has been taken into account. Volume 
fractions investigated in the following are Lp = 0.50, 0.53, 
0.55, 0.57, 0.58, and 0.585. In all states 5 or 10 differ- 
ent systems were prepared from scratch, and 100 or 50 
independent correlation functions were measured adding 
to a total statistics of 500 evaluations of the dynamical 
quantities per state. 

In order to study the specific effects of polydispersity, 
a varying number of bins M has been used to group par- 
ticles according to their size: besides the usual effective 
one-component analysis, il/ = 1, we also discuss a three- 
component, M = 3, and a five-component, M = 5, in- 
terpretation of the data. In these cases, bins of uniform 
width have been chosen. We mainly discuss the collec- 
tive density correlation functions (intermediate scatter- 
ing functions), 

*a/3(9,t) = (6la(<7,i)*(?/3(g)) , (3) 

whose equal-time values are the partial static struc- 
ture factors Sapiq) = '^a[j{q,t = 0). Here, indices 
a,/3 = \^...M label the component bins containing 
^a,/3 particles, and ga{q,t) = N~^^^ Yjk=i exp[«q ■ r^(t)] 
are the collective partial-number-density fluctuations at 
wave vector q, where r^{t) is the position of the k-th 
particle of species a. Note that these correlation func- 
tions are real-valued and depend on q only through its 
scalar invariant q, as the system is isotropic and trans- 
lational invariant. The brackets indicate time-origin or 
canonical averaging in the simulation or theoretical ap- 
proach, respectively. To improve the statistics of this 
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function, a small dispersion in g-modulus was allowed 
for, Sqd = 0.2/^ for qd > 6.0 and Sqd = 0.2/v^ for 
qd < 6.0. The change in Sq can be noted as a kink in 
some wave-vector dependent quantities derived from it 
at qd = 6.0. 

Additionally, we measure the tagged-particle (in- 
coherent) density correlation functions, i.e., the self- 
intermediate scattering function 



(4) 



where g^^{q,t) = exp[iq ■ r^{t)] denotes the Fourier trans- 
formed one-particle density. In this case, averaging over 
all particles of the same species allows to further improve 
statistics in the simulation. {q, t) is connected in the 
\ow-q limit to the mean squared displacement of species 



(5) 



The "polydispersity-averaged" total correlation func- 
tions are recovered by summing over the bins. 



(j){q,t)- - 



S{q) 



(6) 



which in the tagged-particle correlation function reduces 
to 



'(g,t) = (l/M)^$^(g,t) 



(7) 



and the analog expression for the mean-squared displace- 
ment, 



Jr2 = l/A/^5ri(g,i). (8) 



B. Mode-coupling theory 

The mode-coupling theory for M-component mixtures 
builds upon an exact equation of motion derived for the 
matrix of the partial dynamical density correlation func- 
tions $Q,^(g,i) defined in Eq. ([3]). Applying Zwanzig- 
Mori projection operators to the Liouville equation gov- 
erning the microscopic dynamics, one arrives at (sBj 

+ J(g)5-1 ((?)*((?, i) 

t 

+ J{q) I dt' M{q, t - t')^{q, t') ^ . (9) 



with Jajsiq) = q'^SapXav'^Yi a ^^'^ mcmory kernel ma- 
trix Mq{t) which embodies the fluctuating quantities 
and plays the role of a generalized friction coefficient. 
Note that in our case, Vth,a is the same for all pseudo- 
species. Initial conditions for the equations of motion arc 
= 0) =: S{q) and = 0) = 0. We split the 



memory kernel into a contribution describing the regular 
part of the friction, modeled as a Markov process with 
a damping coefficient u = ^ that is chosen in agreement 
with the one taken in the simulations, and a collective 
part describing the slow dynamics. 



M{q,t) « J~\q)u{q)5{t) + M^'^^ {q,t) 



(10) 



where Vajsil) = i^/M with v = 7/2. MCT now approx- 
imates M'^^'^ as a nonlinear functional of the density 
correlation functions, 



M^^^{q,t)^Fmq,t) 
with components 

1 p ^ f d^k 



(11) 



y 

2q XaXp ^-^ 



a'P' 



(27r) 



VaQ'Q"(q, fe,p) 



X <^a'P' (fc, <)*a"/3" [P, t)Vf}p,pn (q, fc, p) . (12) 

The vertices couple the density fluctuations of different 
modes, with q = k + p, 



(qk) 

Vapjiq,k,p) = Cap{k)5a-/ 



(qp) 



q 



q 



Ca^(j))SaP ■ (13) 



Here we have additionally approximated static three- 
point correlation functions in terms of two-point ones. 
Their inclusion is computationally too demanding and 
docs not change the results strongly in dense, non- 
network forming systems such as ours [1^. c{q) is the 
matrix of direct correlation functions defined through the 
Ornstein-Zernike equation 



(14) 



Thus, taking S{q) from simulation, the collective dynam- 
ical density correlators are fully determined in the theory. 

The tagged-particle correlator $^ {q, t) for a singled- 
out particle of species a obeys an equation similar to 
Eq. ®, 

+ / dt'M^^{q,t-t')^iiq,t')=0 (15) 
Jo 

with fl^^q)"^ = q^v^^ ^ and the initial conditions $„((?, t = 

0) = 1, $^((7,i = 0) = 0. We set i/"(g) = v to ob- 
tain damped-Newtonian short-time dynamics. The cor- 
responding memory kernel is evaluated from Eqs. (jl2p 
and (jl3p by considering a (M-l-l)-component mixture in 
the limit of one concentration going to zero: 

M^(g,t)= J-^ [*(<?, t),$'^(g,i)] = 



d^k 



— ^ J (27r)^ 



E 



Vs,a'/3'(9,fc)$a'/J'(fc,0*afe0, (16) 
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with the tagged-particle vertex 

qk 



(17) 



As the memory kernel in Eq. (|16p contains the coher- 
ent correlator, solving Eq. (|15p requires the full coherent 
dynamics to be known. 

By virtue of the expansion = l — q^6r1^^{t)/Q+ 

0(5^), the mean-squared displacement of species a is 
connected to the tagged-particle correlator via 

<Sr2„(i) = lim-|[l-$^(<Z,i)]. (18) 
Its equation of motion follows from Eq. ([T5|) : 



dt5rl^{t) + vt^,^ / dt' AC{t - t')Srl^{t') = &v,,,.,J , 
Jo 

(19) 

with M^{t) = \im.q^oq^M^{q,t) and 



M. 



X $,,;3,(fc,t)$^(fc,t) (20) 



Equation ()20p states that the mean-squared displacement 
is completely determined from the collective and tagged- 
particle density correlation functions; there is no back- 
coupling of the MSD to itself since the phase space at 
q = has vanishing contribution inside the integral. 

Many features of the solutions of the above MCT equa- 
tions are known, especially concerning points asymptoti- 
cally close to MCT glass transitions. We only summarize 
the basic results for completeness, referring to the liter- 
ature H, l4T|-f43j for details. The starting point of the 
asymptotic analysis is to realize that the MCT equations 
allow for bifurcation points for the long-time limit of their 
solutions. Denoting them by 

F{q) = lim *(q,i) F^(g) = lim $'*((?, t), (21) 

one finds that these long-time limits (synonymously 
called glass form factors or nonergodicity factors) are de- 
termined by a set of coupled, implicit nonlinear equa- 
tions, 

Siq)-Fiq)=[S-\q)+T[Fiq)]]-\ (22) 
F'{q) 



F-{q) 



= T^[F{q),F%q)] 



(23) 



In the usual case, it is Eq. that displays bifurca- 
tions: Out of the possibly many solutions to this equa- 
tion, the long-time limit corresponds to the non-negative 
real solution that is largest according to a straightfor- 
ward ordering defined for each q separately and through 
the positive-definiteness relation [44|, . The glass tran- 
sitions of MCT are then the bifurcation points affecting 



this largest solution that arise from smooth changes in 
the control parameters, and the most common case is 
that of a A2 bifurcation where the long-time limit jumps 
discontinuously from the trivial zero solution indicating 
a liquid to a finite value indicating a solid. The solu- 
tion F'' (q) determined as the largest solution of Eq. ((23|) 
will then, in the generic case, inherit the bifurcations of 
F(q), and we will not discuss the possibility of extra sin- 
gularities arising in Eq. ((23)) itself. The generic case is in 
particular obeyed by the problem at hand, quasi-hard- 
sphere tracer particles inside a quasi-hard-sphere system 
composed of particles of roughly equal size. Generally, 
the transition points are then defined as the points where 
the stability matrix C of the nonlinear Eq. ([22]) , given by 



C[H{q)] ■.^2{S^{q)-F^{q))x 
F[F%q),H{q)] {S%q)-F^{q)), 



(24) 



has a unique critical eigenvector H[q) with eigenvalue 
unity. H{q) is also called the critical amplitude (up to 
trivial normalization factors that are sometimes split off 
from it). We will denote quantities corresponding to such 
a transition point with superscript, e.g., F^ and -Fg''^. 

On the liquid side of the transition, correlation func- 
tions follow a two step relaxation scenario: for times large 
compared to the characteristic time of the short-time mo- 
tion, t ^ to, they decay with a time fractal ~ to the 
plateau, which extends until the /3-relaxation time scale 
to-- For t ^ tcr, the decay from the plateau sets in with 
the von Schweidler law, ^ —t^, initiating the final a re- 
laxation that is characterized by a second time scale t'^. 
The asymptotic analysis proceeds by analysing the equa- 
tion of structural relaxation, where time-derivatives that 
affect only the short-time motion have been dropped, 

*(9,t) = S{q)M{q,t)S{q) 

- j^S{q) j\t'M{q,t-t')^{q,t'). (25) 

Identifying the distance of the correlator to its plateau 
value as a small parameter a, one extracts the two time 
scales that diverge upon letting cr — > 0, 



(26) 



where 7 = l/(2a) + 1/(26), and a > and b > are 
nontrivial and nonuniversal exponents determined by the 
details of the interaction potential (see below). The sep- 
aration parameter a is, in leading order, linearly con- 
nected to the change in control parameters, viz. in our 
case cr = Ce + 0{e^) with e = {(p — ip'^^j^f . By conven- 
tion, e < indicates a liquid state, e > the glass. Note 
that t'^jtcr also diverges as 0, so that asymptotically 
close to the MCT transition, an ever larger window for 
structural relaxation around the plateau open. In prac- 
tice, this window is cut short for large — e by preasymp- 
totic corrections, and for small |e| when the theory fails 
to describe residual ergodicity-restoring processes in the 
glass. 
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For times t = t/t„ where the correlator t) is close 
to F^{q) one makes the following Ansatz 



*(<7, K) = F\q) + ^G{q, i) + 0{\a\) , 



(27) 



and the uniqueness of the critical eigenvector at the bifur- 
cation point implies the factorization theorem, G{q, t) = 
Hq ■ G{i). The function G{t) is determined by the so- 
called /3-scaling equation 



f dt'G{i - t')G{t') - A(G(£))2 + sgna = . 
dt Jo 



(28) 



The nonuniversal details of the vertices enter in this equa- 
tion only thorugh the exponent parameter. 

A = Hq:{S%q)-F-{q))J^%q)[H{q),H{q)]x 

(S^q) - F^q)) /N (29) 

where Af = H{q) : (ff(g)(S^(g) - F%q)r^Hiq)), 
and the double-dot operator includes contraction over 
q. Here, H(q) is the left-eigenvector corresponding to 
H{q). 

For times <^ t <^ the decay to the plateau at 
the critical point is then governed by the /3-relaxation; in 
leading order, 

*(g, t) = F^iq) + H{q){tlt,)-^ + O ((t/t.)-^'^) (30) 

where a is determined as solution 0<a; = a<l/2of 

^-r(i-2x)- ^^^^ 

For times ^ t to-' and for cr < the decay of the 
correlator is described by the von Schweidlcr law 

*(q, i) = F^{q) ~ H{q){t/t„f + O {{t/U)'') . (32) 

Here b, the von Schweidler exponent is determined from 
the negative solution 1 > 6 = —x > of Eq. ([3T|) . 

In the glass (cr > 0) the nonergodicity parameters be- 
have like 



F{q)=F^{q)+H{q) 



1- A 



0{a). 



(33) 



Again we define the "polydispersity-avcraged" total non- 
ergodicity parameters and the total critical amplitudes 
by summing over the bins, 

For timescales t = t/t'^ ~ 1 and f ^ Lp'^ the a-master 
equation can be derived 

Mq,i) = S{q)M{q,i)S{q) 

- 4^(g) / dt'MiqJ^t')4>{q,t'). (35) 
dt Jo 

This equation states that all correlators should collapse 
on the same function when rescaled by an appropriate 
scaling time i. This is due to the invariance of Eq. ([55)1 
when rescaled in time and is the mathematical manifes- 
tation of the time-temperature superposition principle. 
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FIG. 1: (Color online) Total static structure factor extracted 
from the simulation. The packing fractions ip = 0.45, 0.5, 
0.55, and 0.57 are marked with crosses (red), squares (green), 
circles (blue) and triangles (magenta), respectively. The black 
solid line is the Percus-Yevick result for tp = 0.516. The inset 
shows a magnification of the peak. 



Numerical details of the MCT solution 



MCT calculations were performed with static struc- 
ture factor matrices taken from the simulations; to access 
packing fractions for which no simulations have been run, 
linear interpolation in ip was used for S{q). As in the 
simulation, calculations with a one-component (M ~ I), 
and two different multi-component binnings. A/ = 3 
and M = 5, have been performed (see Sec. Ill A|) . To 
choose effective diameters da for the different bins, we 
have followed the moment approximation jssl [36j : for 
M = 3, the three values {di, (i2, da} and equal concen- 
trations allow to match the actual polydispersity dis- 
tribution up to the second moment (with the require- 
ment that the skewness of the distribution vanishes), 
and for M = 5, additionally the fourth central mo- 
ment can be matched using equal concentrations. This 
leads to the choice {di, (i2, ds} = {1 — l/\/200, 1, 1 + 
l/\/200} in the M = 3 system, and {di, d2, da, d4, ds} = 
{0.91675, 0.96254, 1, 1.03745, 1.08325} for the M = 5 sys- 
tem. 



The problem of solving Eq. © in full has to be tack- 
led numerically, choosing a suitable discretization. The 
wave number grid selected here is g € [0.1, 60.0] in steps 
of Aq = 0.2. The discretization is a compromise be- 
tween calculation time and being as close as possible to 
the simulation structure factors. The initial time step 
was chosen as 5t = 10~^. After every 128 steps the step- 
size was doubled in order to cover logarithmically large 
intervals in t. 
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M = 1 0.566 
M = 3 0.537 
M = 5 0.535 
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FIG. 2: (Color online) Partial static structure factors ex- 
tracted from the simulation by binning the polydisperse sys- 
tem into M = 3 species, at packing fraction (p = 0.5. The 
upper curves show the diagonal terms Sii{q) (red crosses), 
S22{q) (green squares), 5*33 (g) (blue circles). The lower curves 
show the off-diagonal terms Si2{q) (red crosses), Si3{q) (green 
squares) and 5*23(9) (blue circles). 



III. DATA ANALYSIS 
A. Structure factors 

Wc first address the static structure factor that serves 
as an input to MCT, S{q) and S{q) — J2afS ^apiq) ob- 
tained from the simulations. In Fig. [U the averaged 
structure factor S{q) is shown for different packing frac- 
tions. It exhibits the standard features known for dense 
liquids where excluded volume is the main interaction 
effect: a pronounced first peak is visible that shifts to 
higher g- values and increases in intensity with increasing 
packing fraction, reflecting a decreasing average next- 
neighbor distance and increased local ordering in the 
denser system. The figure also includes S{q) as calcu- 
lated from the Percus-Yevick (PY) approximation for 
monodisperse hard spheres [1^. It is known that this 
approximation performs quite well, but somewhat over- 
estimates the amount of ordering present in the system; 
this is expressed by a shift in density values. Still, com- 
paring densities where the first main peak is reasonably 
well described within PY, the amplitude of the second 
peak in S{q) is notably overestimated by the approxima- 
tion in comparison to our simulation results. Note that 
in MCT, the main contribution to the memory kernel 
causing slowing down and arrest comes from these am- 
plitudes in S{q) — 1. As larger wave numbers contribute 
more strongly to the three-dimensional integral, not just 
the main peak of S{q), but also the shape of its large-q 
envelope determine the MCT dynamics; one can thus an- 
ticipate from Fig. [T] that MCT calculations based on the 
PY structure factor will clearly underestimate the critical 
packing fraction. 

In Fig. [5J we show the partial structure factors result- 
ing from a binning of the simulation data into a three- 



TABLE I: Table of MCT-calculated critical packing fractions 
for the different polydispersity models using M bins. Cal- 
culations are based on the computer-simulated partial static 
structure factors. 



component system. All diagonal terms as well as the off- 
diagonal terms (except for a trivial shift by the constant 
l/(xQ,a;^)^/^ = 1/3) arc very similar as expected for a 
mixture of almost equal constituents. The slightly differ- 
ent average next-neighbor distances for the different par- 
ticle sizes cause corresponding shifts in the oscillation fre- 
quencies and hence the positions of the peaks in Sap{q)- 
This effect is most pronounced in the diagonal terms, as 
these only contain information from one distinct particle- 
size bin. Adding these partial structure factors recovers 
the S{q) shown in Fig. [T] and elucidates that the com- 
paratively weak high-q peaks visible there arc the result 
of a destructive interference of the slightly shifted oscil- 
lations in the multi-component S{q). However, the main 
peak is strong enough to be less affected. Thus, for sys- 
tems whose interactions are close to hard-core repulsion, 
such as our system, averaging S{q) mainly results in an 
underestimation of the coupling strength, i.e., too fast 
dynamics in the one-component calculation as compared 
to a multi-component calculation. Indeed, our MCT cal- 
culations discussed below essentially confirm this expec- 
tation, that in the range of q around the first peak in S {q) 
and above, the changes between pre- and post-averaging 
are minor once the shift of the critical packing fraction 
is acocunted for. Yet, for small the situation is more 
intricate. Let us also note that for systems whose glass 
transition features are dominated by the large-q behavior 
of the static structure factor, the destructive interference 
induced by pre-averaging the MCT input can lead to se- 
vere changes in the qualitative dynamics (2l| . This par- 
ticularly concerns polydisperse colloid-polymer systems 
with short-ranged (depletion-induced) interaction among 
the colloids. 



B. Critical point density 

The critical packing fraction Lp'^ of the MCT glass tran- 
sition was calculated with a bisection algorithm which re- 
turns the point where the eigenvalue of the stability ma- 
trix C , Eq. (p4)) , is unity. Table U summarizes the values 
of (^^ along with the corresponding exponent parameters 
A and the MCT exponents corresponding to these values, 
for the M = 1, 3, and 5 calculations we performed. 

As expected from the discussion of the static structure 
factor above, the glass transition point in the calcula- 
tions shifts to lower densities with increasing M. Re- 
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markably, the most pronounced change aheady occurs 
between M = 1 and M = 3, whereas the further increase 
to AI ~ 5 docs not alter ip'^ or A substantially. We thus 
consider M = 5 to be already close to the limit of many 
components that should in fact be taken to describe a 
truly polydisperse system. At this point we would like 
to stress that the downshift of the glass transition point 
for M = 1,3,5 is only related to the different MCT- 
models used here, as the polydispcrsity in the simulation 
is not changed. It should be noted that still, even the 
value V3^^5 ~ 0.535 is significantly above the value cal- 
culated within the PY approximation, ippy « 0.516 [4^ . 
However, the MCT glass transition from an asymptotic 
analysis of the simulation data yields a value ifMD 

« 0.59 

[2^, much higher than the values calculated within the 
theory even with the correct static structure information, 
but in reasonable agreement with the commonly quoted 
value for polydisperse hard-sphere like colloidal suspen- 
sions, (Pcxp ~ 0.58 [IJ. This is a well known shortcom- 
ing of MCT, and we consider the better agreement we 
obtain for the one-component calculation, V'm=i ^ 0.566 
as fortuitous. Note also that one usually expects poly- 
dispcrsity to shift the glass transition to higher densities, 
as the overall packing efficiency incresases |48| . However, 
the role played by the different shapes of the polydispcr- 
sity distributions is not well studied. Within MCT, this 
commonly observed trend of increasing ip"^ in mixtures is, 
at least for binary hard-sphere mixtures, only predicted 
for size ratios S < 0.7, while for S closer to unity, the 
inverse trend is found [36|. Relating S to the extreme da 
in our calculations, even our M = 5 system only corre- 
sponds to 5 ss 0.85. 



C. a-process analysis 

Next, we test the validity of a. scaling for our simula- 
tion data. According to Eq. ([55)1 . plotting correlators as 
functions of tjt'^ should collapse the data for long times, 
with a master curve extending s& ip ip'^ from below. 

To determine a relaxation time r cx t'^ from the data 
alone, we have shifted the correlation functions at a single 
fixed value oi q~ T.i/d to coincide with the correspond- 
ing curve at the highest packing fraction in the liquid, 
(f = 0.585, at long times. After this procedure, the va- 
lidity of the scaling can be checked by requiring that r 
be independent on q. Note that qd « 7.3 corresponds 
to the position of the main peak in the averaged static 
structure factor. We have chosen this value since here 
the strength of the a process is maximal, and the best 
separation from the (3 regime is achieved. Fig. [3] shows 
the result of this scaling for four different wave numbers. 
An a-master curve clearly is approached, with the scal- 
ing regime for the highest two densities extending over 
about two orders of magnitude in time. The strong cou- 
pling of the a-relaxation on local length scales predicted 
by MCT is observed, as the same scaling factor r rescales 
the correlators for different wave vectors. 
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FIG. 3: (Color online) Simulation correlators rescaled by the 
a timescale t oc t'^ for four different wave numbers q as in- 
dicated, where r was determined by shifting the results for 
lower ip to agree at long times with the ip = 0.585 curve at 
qd = 7, and does not depend on q. Packing fractions shown 
are ip = 0.585 (black plus symbols), ip = 0.58 (red crosses), 
p = 0.57 (green stars), ip = 0.55 (blue open squares), p = 0.53 
(magenta filled squares) and p = 0.50 (black circles). 

In the lower panel of Fig. |3l the a scaling is exhibited 
for small qd; here, deviations are visible at qd ~ 3, where 
the simulation data does not exhibit a common shape for 
the a-relaxation part of the correlators. The highest p 
shown indicate a decay that is either more stretched or 
exhibits a further infiection point in the 0(t)-versus-logi 
plot below the plateau. We will return to such features 
below in the discussion of hydrodynamic intcrdiffusion 
modes that exist in multi-component systems and inter- 
play with the structural relaxation at low qd. 

A common description of the shape of the a relaxation 
is in terms of stretched-exponential (Kohlrausch) laws, 

^q,t)^Agexp[-it/Tgf-]. (36) 

Here, /3q is the stretching index, required to be /3q < 1 
for structural relaxation in equilibrium systems. While 
the a-master function from MCT is different from the 
Kohlrausch form, the theory predicts that for large wave 
numbers, the two functional forms become identical, and 
/3g_j.oo — >■ b [i^. Tq is commonly referred to as the a- 
relaxation time, and it is connected by a g-dependent 
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FIG. 4: (Color online) Fits of stretched-exponential 
Kohlrausch functions (dashed black lines) to the simulated 
coherent density correlators (circles) at (p — 0.585. The q- 
values are from top to bottom qd = 6.6 (red), 7.4 (black), 9.8 
(blue), 12.8 (green) and 15.6 (magenta). The fit range was 
chosen as t £ [100 : 10^]. Solid black lines show von Schwei- 
dler fits up to second order, Eq. (|37p . with = 1000 and 
6 = 0.53, fitted in the range t € [20 : 336]. 



but density-independent prefactor to the scaling time t'^ 
appearing in MCT. The Kohlrausch amplitude Aq < 1 
can be taken as an estimate of the MCT plateau value 
fq, and since the a process by definition starts below this 
plateau, Aq < fq is expected. In practice, however, the 
separation of the a process from the /? relaxation is often 
not clear enough to warrant this restriction. 

In general, Kohlrausch fits arc hindered by some sub- 
tle problems that are often overlooked. Lacking a clear 
separation of the a process, the fit parameters exhibit a 
dependence on the fit range. It is not clear a priori how 
to choose the optimal fit range, as for very long times, 
one expects the relaxation to become (non-stretched) ex- 
ponential again (and to be covered within the noise of 
any experiment), and for short times, deviations are seen 
that can be understood within MCT to be due to the 
difference between /3q and b for finite q. 

We have tried to fix the fit range of our stretched- 
exponential fits such that the fit parameters exhibit only 
a weak dependence on the fit boundaries. Examining 
the q = Qp correlators for (p = 0.585, this leads to 
t £ [10^, 10^]. In Fig. m examples of such fits are shown. 
Although the agreement is generally convincing, let us 
stress that these fits are a pure data analysis, not taking 
into account the full numeric MCT calculations we will 
discuss below. The curves in Fig. |3] show what one can ex- 
tract parameters like the plateau value or the von Schwei- 
dler exponent b without having performed full MCT cal- 
culations. In Fig. 131 we additionally show the results 
of von Schweidler fits to the simulation data: adapting 
Eg. and extending it by the next-to-leading order 
[43, we fitted 

$(q, t)^fq- hq{tlt„f ■ (1 - kq{t/t,f) , (37) 
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FIG. 5: (Color online) MCT results for the critical noner- 
godicity parameters obtained by binning the simulated 
static structure-factor data into A4 — 1 (black dashed), 
A'l — 3 (red solid) and M = 5 components (blue dash-dotted). 
The grey dotted curve is the MCT solution with the Percus 
Yevick static structure factor (M = 1). Circles show the 
plateau values of the simulation curves obtained by fitting 
Kohlrausch functions to the coherent simulation correlators 
at 1^ = 0.585 (Fig. 2]). Triangles are obtained by von Schwei- 
dler fits (Eq. (|37l) and Sec. [HID]). 



determining fq, hq, and kq by fitting and fixing b = 0.53 
to be consistent with the analysis of a-relaxation stretch- 
ing presented below (Fig. [5]). Note that a free fit of 
von Schweidler 's law without fixing b is usually ambigu- 
ous as the determination of the MCT exponents directly 
from data bears uncertainties [HO]. Eq. ([57]) has been 
used in the range t S [20; 336]. The von Schweidler re- 
sults then represent about two decades in time of the 
correlation functions. 

Figure [5] shows the results for the amplitude Aq, 
together with the nonergodicity parameters fq calcu- 
lated from MCT, and the estimates of fq obtained from 
von Schweidler fits. The values obtained from all three 
methods in general show good agreement, although some 
details warrant discussion. Let us first turn to the large-g 
regime, qd > 7. Here, the von Schweidler fits, Eq. (|37p . 
yield fq that are in good agreement with the Kohlrausch 
amplitudes Aq, and the relation Aq < f^ is reasonably 
well fulfilled. The MCT calculations somewhat underes- 
timate fq, although the agreement is soupcon improved 
for the multi-component systems with M = 3 and M = 5, 
in particular for values of q where the averaged static 
structure factor (and hence also the /^-versus-g curve) 
shows minima. The underestimation of fq by MCT can 
be attributed to the underestimation of ip'^: MCT de- 
scribes arrest at lower densities, but fq increases with 
density as the denser glass is stiffer with respect to den- 
sity fluctuations. For comparison, we also show in Fig. [5] 
the fq obtained by employing the FY approximation 
within MCT; this result agrees well with the A/ — 3 
and M = 5 calculations, indicating in particular that for 
the arrest of small-wavelength fluctuations, polydisper- 
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FIG. 6: (Color online) MCT results for the tagged-particle 
critical nonergodicity parameters f^''^ for the M = 1 (dashed), 
M = 3, (solid) and M = 5 (dash dotted) polydispersity mod- 
els. The solutions for M = 3 and M — 5 cannot be distin- 
guished on the scale fo the figure. Open circles and triangles 
show plateau values determined from fitting Kohlrausch and 
von Schweidler functions, respectively, as in Fig. [S] 



sity and particle softness do not play a major role. One 
in fact recognizes that the truly monodisperse calculation 
(with PY input) is in better agreement with the polydis- 
perse data than the one-component calculation using the 
polydispersity-averaged S{q). This can be intuitively in- 
terpreted: pre-averaging in S{q) artifically reduces large- 
q static correlations, but fq is in essence a dynamical 
quantity, and the effect of the dynamics on these large-q 
correlations needs to be included to describe the degree 
of dynamical arrest. Still, in Fig. [5] this is a merely quan- 
titative effect. 

For the large-wavelength regime, qd < 6, the situa- 
tion is more differentiated: here, both one-component 
calculations differ more notably from the results obtained 
by pure data fitting, and the results calculated with the 
M = 3 or M = 5 theory. In particular, the PY-based 
one-component result shows a rather weak g-dependence 
of the plateau value for small q, and yields « 0.4 in this 
regime. In contrast, the multi-component results show an 
upturn of as q is decreased, with plateau values rising 
to almost 0.8 at the lowest q accessible in the simula- 
tion. The one-component MCT calculation based on the 
simulated (prc-averagcd) structure factor also shows this 
upturn, but less pronounced. 

This behavior can be rationalized by the effect of 
polydispersity: first, going over from the strictly one- 
component PY approximation to the simulated S{q) al- 
ready contains static, pre-avcraged information about the 
size distribution, and this is sufficient to explain qualita- 
tively that polydispersity tends to increase the stiffness 
of the frozen structure towards long-range density fluctu- 
ations. The post-averaged calculation is needed to cap- 
ture this trend quantitatively. It shows that additionally 
the freezing out of the interdiffusion process increases the 
average fq at small q (5l| . 
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FIG. 7: (Color online) a-relaxation times r, obtained by 
stretched-exponential fits to the MCT master curves for the 
M = 1 (black dashed), M = 3 (red solid), and M = 5 (blue 
dash-dotted) polydispersity approximations. All values have 
been rescaled by their value at qpd = 7.3, r* — Tq/rq^. The 
upper panel shows data corresponding to the collective den- 
sity correlation functions, the lower panel those corresponding 
to the tagged-particle analog. Circles correspond to the Tq ex- 
tracted from the simulation data at tp = 0.585, as in Fig. [31 



The tagged-particle analog of the quantities shown in 
Fig- El fq''^i is exhibited in Fig.|6l Here, the differences in 
the different methods of determining the plateau height 
are much smaller, as intuitively expected; the incoher- 
ent /*''^(g)-versus-g shapes do not show the oscillations 
typical for the collective quantities. Only at the largest 
wave numbers investigated here, qd > 12, some differ- 
ences between the MCT calculation and the fit results 
become apparent. Here, the theory again underestimates 
the amount of arrested density fluctuations, and this can 
agian be rationalized by its lower critical density. For 
small q, all f'^''^{q) have to approach unity, so that any 
differences are trivially wiped out. 

Let us now turn to a discussion of the g-dependence 
of the relaxation time Tq and the stretching Pq resulting 
from Kohlrausch fits. To obtain equivalent values also for 
the MCT curves, we have fitted the theory's a-master 
curves with stretched exponential functions. Numeri- 
cally, the master curves have been approximated by cor- 
relators close to the transition point: choosing e — — 10~^ 
in Eq. (j9l) proves sufficient to effectively solve Eq. ((25)) 
for the timescales of interest. For the simulation, we 
again restrict the discussion to the highest density avail- 
able. If = 0.585, which however corresponds to a different 
e. Thus, to enable a meaningful comparison, relaxation 
times are scaled to coincide at q^d ~ 7.3. 

Figure [7] shows the resulting a-relaxation times for 
both the coherent (upper panel) and incoherent (lower 
panel) dynamics. A similar distinction into two regimes 
as above arises: for qd > 6, the g-dependence of the re- 
laxation time is excellently predicted by the theory. For 
Qd < 6, this holds for the collective correlation functions 
only when comparing with the multi-component calcu- 
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FIG. 8: (Color online) Kohlrausch stretching parameters Pq 
as functions of wave number q, determined from fits to the 
MCT Q-master curve for the M = 1 (dashed), M = 3 (solid), 
and M — 5 (dash-dotted) polydispersity moment approxima- 
tions. Pq determined from fits to the simulation data (as in 
Fig. El, are shown as open circles. The upper (lower) panel 
shows results from analyzing the collective (tagged-particle) 
quantities. 



lation. Even the pre-avcraged one-component calcula- 
tion is in qualitative difference, since there, Tq^o ap- 
proaches a constant, resulting in a weak g-dependence 
for all qd < 6. Instead, in the multi-component theory, 
Tq ~ I/q^ as g — > for every single matrix element of 
the partial-density correlation-function matrix. This has 
a clear physical interpretation: Tq ^ const, reflects the 
fact that the overall momentum of the system is con- 
served. However, this is not the case if one considers a 
single species inside a mixture of components only, since 
momentum can be exchanged between the species. Only 
the sum of the partial momenta is thus conserved, which 
is reflected in the existence of an appropriate zero eigen- 
vector in the MCT memory kernel; however, this does not 
transcede to the averaged correlator defined by Eq. © 
itself. Strikingly, the resulting rise in Tg as g — >■ that 
is apparent in the M > 3 MCT calculations is in very 
good agreement with the simulation results, with only 
the lowest q values of the simulation deviating slightly 
(which may be due to minute finite size effects). 

For the same reason (momentum conservation or 
rather the lack thereof for a single species), ~ l/q^ 
is expected in any system. Indeed, this is seen in all 
curves in the lower panel of Fig. [T] However, the diver- 
gence predicted by MCT is stronger than the one seen in 
the simulation, an effect visible even at finite q. Only in 
the large-g regime defined above is the MCT description 
of the t| quantitatively accuracte. At the largest q ex- 
tracted from the simulation, the theory in turn somewhat 
underestimates the relaxation times. 

Turning to the stretching exponents Fig. |8l the 
agreement between MCT and simulation is less favor- 
able. The values obtained from fitting the theory curves 
are systematically too high, an effect that will become 



evident also below when discussing the full correlators. 
Only the qualitative behavior of Pq with q is qualitatively 
in agreement, although the difficulty of determining /3g in 
the simulation for large q, where the amplitudes Aq are 
already rather low, does not allow a detailed discussion. 

Since the a master curve strictly is a Kohlrausch func- 
tion only in the limit g — >■ oo [49t , emphasis should mainly 
be placed on the behavior of (3q at large q, and at small 
q, as we will discuss below. For large q, the fits to 
the MCT curves nicely exhibit the asymptotic behav- 
ior, Pq^ao b, with a value of 6 sa 0.6, consistent with 
the values given in Table HI Also the simulation-fitted 
stretching exponents are compatible with the approach 
to a finite constant at large q, albeit with a somewhat 
lower value, b « 0.5; we can take this difference as an 
indication for the error inherent in the value of A as cal- 
culated within MCT, in particular since the value of b is 
in good agreement with the b ~ 0.53 that results from 
an independent /3-relaxation analysis of the simulation 
data (see below). For the simulated incoherent corre- 
lation functions, however, no such large-g limit can be 
identified for (3^; the reason for this behavior is unclear. 

At low g, the appearance of a diffusion mode in the in- 
coherent correlator demands — >■ 1 for g — >■ (since hy- 
drodynamic relaxation functions are just exponentials). 
The fits to both the different theory calculations and to 
the simulation data confirm this. For the coherent /3q, no 
such statement holds in the true monodisperse system; 
there is no collective single-component diffusion mode. 
Hence, fSq^o < 1 is found for the MCT calculation based 
on the pre-averaged structure factor. For the same rea- 
son discussed in connection with the Tq ~ 1/g^ behavior 
above, however, the /3g from fits to the simulation data 
do exhibit an increase as q decreases towards zero. In 
principle, a similar trend should be found in the multi- 
component MCT calculations. However, a peculiarity 
arises here, that prevents us from determining meaning- 
ful values of /3g in this case for qd < 4. Here, the appear- 
ance of a number of distinct interdiffusion processes typ- 
ical for a multicomponent system leads to a-relaxation 
curves that arc superpositions of a small number of ex- 
ponentials, leading to master curves that exhibit multi- 
ple "shoulders" ; for the corresponding relaxation spectra, 
this corresponds to multiple a peaks [sii . i52i] . Such in- 
terdiffusion processes in principle also exist in the simu- 
lation; however in this truly polydisperse system, a large 
number of them is combined to a relaxation curve that 
is again reminiscent of a single a process (akin to a het- 
erogeneous superposition of single relaxators). 



D. /3- process analysis 

We complete the asymptotic analysis of our simulation 
data by investigating the /3-scaling regime. On approach- 
ing ip'^, MCT states that the correlation functions be de- 
scribed in leading order by Eq. (|27p . However, testing 
this relation involves a number of fit parameters whose 
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FIG. 9: (Color online) /? analysis of the simulation data at 
— 0.585. Functions X{q,t) calculated from Eq. (|38|l by 
fixing t' = 10.0085 and t" = 52.5810 are shown for the collec- 
tive (tagged-particle) correlators in the upper (lower) panel. 
Different wave numbers q were chosen as labeled. The MCT 
factorization theorem is validated by observing data collapse 
for different q in a time window spanning [t' ,t"], and by the 
/3-master curve shown as a dashed line. 



determination is difficult. An approach to testing the 
first main prediction of MCT for the /3-relaxation win- 
dow, viz. the factorization theorem, is to consider the 
function [SOl 



X{q,t) 



<l>{q,t) - ^{q,t') 
4>{q,t') - ct>{q,t") 



(38) 



with times t' and t" fixed to be inside the scaling regime. 
Equation (^T]) then predicts X{q^t) = xiG{t) — X2 to be 
independent on wave number; hence, superimposing the 
functions X{q, t) for different g, one should be able to fix 
the two times t' and t" uniquely such that a time window 
appears where all X{q^ t) collapse. The procedure has the 
advantage that the critical amplitude drops out and thus 
does not need to be determined by fitting. It was shown 
to work very reliably in a binary Lennard- Jones mixture 
by Gleim and Kob |3G |. 

As shown in Fig. [HI such a collapse is indeed possible for 
the simulation data ai ip ~ 0.585, where we have chosen 
t' = 10.0085 and t" = 52.5810. Both the collective and 
the tagged-particle density correlators collapse for all the 
wavenumbers investigated for a region spanning [t' ,t"] 
and slightly extending to both smaller and larger times. 
Estimating that e = -0.015 (-0.017, -0.020) for the 
M = 1 (M = 3, M = 5) analysis, one cannot expect the 
first-order asymptotic result for the /3-relaxation function 
to hold over more than one decade in time p^ : indeed 
this is roughly what we observe. 

A stronger test of the MCT asymptotics implicit in 
Fig. iniis the so-called ordering rule: since in the next-to- 
leading order corrections to the factorization theorem the 
same g-dependcnt correction amplitudes appear both for 
the early-time deviations and for the long-time correc- 
tions , correlators that lie, say, above the /3 correlator 



for short times must also deviate in that direction for 
long times. Thus, numbering the correlators in the or- 
der in which they deviate from the asymptote for short 
times, the same numbering should be found on the long- 
time side. Figure |9] confirms this. As was also found in 
Ref. Il^l for the MCT calculations based on the FY struc- 
ture factor for hard spheres, this ordering rule is even 
preserved among the correlators at long times, when the 
/? correlator already violates it. This effect can be seen in 
Fig.iniby noting that the leading-order asymptote (drawn 
as a dashed line) emerges between the qd = 12.8 and 
qd = 15.6 correlators, but already crosses the curves for 
smaller qd at around t sa 500, while the correlators obey 
the ordering rule up to the time window plotted, t < 10^. 

In order to extract the critical amplitude hq from 
the simulation, one can define a function in analogy to 
Eq. dSll) by 



(l){q,ti) - (I){q,t2) 



(39) 



with ti, t2 chosen in the /3-scaling regime. The last equal- 
ity follows again from Eq. ([27|) and thus allows us to ex- 
tract the critical amplitudes up to a factor hq^. Since 
Eq. ([39]) becomes independent on the times chosen as 
long as they are in the /3-relaxation window, we can fur- 
ther improve the statistics of Yq by averaging over two 
time windows 15 3i 



E"=?'/'(9,ii) -E"=„/2+i<^('?'*j) _ hq 



Y„ 



Ei=l -^(^O, tj) - Ej=„/2+l <^(90, tj) ^' 



(40) 



'90 



and use all the data points tj within the /3-scaling 
regime, which leads in our case for — 0.585 to tj £ 
[10.0085; 52.5810]. 

As a cross-check, we have also determined the critical 
amplitudes by directly fitting the von Schweidler expres- 
sion including its leading-order correction to the late /3 
regime, as described in conjunction with Eq. p7p. Fig- 
ure [TO] shows the results for the critical amplitude of 
the coherent density correlators, hq. Reassuringly, both 
determinations of hq give results that are fully consis- 
tent with each other. Also shown in the figure arc the 
MCT-calculated amplitudes. For all three values of the 
number of components chosen, Af , the data are in very 
good agreement, with strongest deviations setting in for 
qd < 6. In particular, the strong dip in hq around q~qp 
is well reproduced. In general, the shape of the /i-versus-g 
curve is in this regime quantitatively captured alread y by 
the Pcrcus-Yevick approximation discussed in Ref. [43 |. 
At small wave numbers, qd < 6, deviations set in that are 
the analog of those discussed above in connection with 
and Tq: polydispersity affects the long- wavelength 
fluctuations in the system, and using the pre-averaged 
static structure factor within MCT cannot describe these 
mixture-specific features. It is intuitively clear that, since 
the actual is larger than its one-component estimate 
at small qd, the opposite has to hold for hq, as the nor- 
malization of the correlation function implies fq + hq < 1. 
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FIG. 10: (Color online) Critical amplitudes hq calculated 
within MCT for the M = 1 (dashed), M = 3 (solid), and 
M — 5 (dashed-dotted) moment approximations to the sim- 
ulated polydispersity distribution. For comparison, the PY- 
based theoretical result is also shown (cyan line). Open cir- 
cles and triangles mark the amplitudes determined from the 
ifi = 0.585 simulation data via the function Y{q), Eq. (|40p . 
and via the von Schweidler fits discussed in conjunction with 
Fig. [SI respectively. The results for Y{q) have been scaled 
by a factor 0.61 to account for the unknown normalization in 
this procedure. 




FIG. 11: (Color online) Critical amplitudes hq calculated 
within MCT for M = 1 (dashed), M = 3 (solid), and M = 5 
component (dashed-dotted) approximations to the polydis- 
persity distribution, and by using the Percus-Yevick static 
structure factor (cyan). Open circles are the corresponding 
amplitudes determined from Y{q), Eq. (|40p . scaled by a fac- 
tor 0.63. Triangles show results from von Schweidler fits as in 
Fig. [5] 

In Fig. [11] the tagged-particle critical amplitudes /i^ 
are displayed. Again, the two procedures to determine 
this quantity from the simulation data alone agree. The 
result shows a peak around qd « 10, while for g — >■ 0, 
/i* — >■ follows from hydrodynamic laws. In contrast to 
the coherent amplitude hq, for the MCT results show 
more pronounced deviations from the simulation values. 
The theoretical quantities exhibit only a weak depen- 
dence on the number of component bins Af in this case. 



and peak around qd « 12.5, i.e., at slightly larger wave 
numbers than what is observed in the simulation. Gen- 
erally, a shift of the /i^-versus-g curve arises, indicating 
that MCT gives a wrong estimate of the relevant length 
scale for the tagged-particle motion in the (3 regime. It is 
unclear whether this mismatch can be attributed to the 
mismatch in critical packing fractions ip"^ between theory 
and simulation, as was done for f^''^. It is also notable 
that the disagreement is pronounced only in the tagged- 
particle critical amplitude; absorbing it into an effective 
wave number, as done in Ref. [29j . would in fact worsen 
the agreement for the collective amplitude as can be seen 
in Fig. [lOl 

IV. FULL MCT-ANALYSIS 

Having established the generic MCT scenario for the 
simulation, we now present the numerical solutions of 
the full (non-asymptotic) MCT equations. In principle, 
the dynamical correlation functions thus obtained can 
be directly compared to the corresponding quantities ex- 
tracted from the simulation. However, the mismatch in 
the values neccessitates a comparison not at equal 
densities, but at, in principle, equal separation from the 
respective transition points. In the spirit of the MCT 
asymptotics, a comparison should involve matching the 
separation parameter cr; however, this is not easily de- 
termined for the simulation, since the true functional de- 
pendence between a and the control parameters is not 
known. Only asymptotically close to the transition do 
we have a oc e, with a pre-factor that also can only be 
calculated within MCT (and thus might be in error). It 
is therefore practical to perform a fitting of the pack- 
ing fractions used in MCT, ipucT to the nominal ones 
used in the simulation, (p, for each of the MCT systems 
with a different number of components M. We have per- 
formed this fitting based only on the coherent correlators 
at q = qp] the comparison for all other wave numbers, and 
for all tagged-particle quantities then is parameter-free. 

It should be noted that a similar fitting procedure was 
already performed in Ref. [l^ for the tagged-particle 
data alone. There, however, it was found that an error in 
the relevant length scale (as discussed in connection with 
Fig. [TT|) could be absorbed by adjusting also the values 
of q in the comparison. Such a procedure effectively al- 
lows to improve the agreement for the plateau values in 
the incoherent correlators, since the /"''^{q) are monoton- 
ically decreasing with increasing q. In the present case, 
no such shifting of wave numbers is allowed for, since for 
the collective f^iq), no such argument holds. 

A. Collective Dynamics 

MCT ascribes the dramatic slowing down in the col- 
lective dynamics approaching a glass transition to a bi- 
furcation scenario, where upon smooth variations of all 
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FIG. 12; (Color online) Packing fractions vpmct chosen for the 
five-component (triangles), three-component (circles), and 
one-component (squares) MCT fits to the simulation corre- 
lators, as functions of the simulated packing fraction ip. The 
dashed lines are linear regression fits, (^mct = cup 4- 6, with 
parameters a = 1.0001 (0.9506, 0.9083) and h = -0.0285 
(0.0288, 0.0288) for M = 1 (M = 3, M = 5). Crosses 
mark <^mct from a fit to the mean-squared displacement only 
(M = 3). Horizontal dotted lines show the MCT predictions 
for the critical point for each of the three models; their inter- 
sections with the linear regression curve marks the estimated 
p'^ for the simulation, as noted by arrows on the abscissa. 



control parameters, a qualitative change in the solutions 
occurs at long times. In fitting the control parameter 
Vmct, it is therefore essential to check that the relation 
^^MCT i}P) does not show signs of singular variation itself. 
In the ideal case, this relation should be linear, as long 
as one considers density intervals where is still large 
compared to the shift A(p'^ = i^ipf — ^mct)- Then, such 
a relation indicates that the a values calculated from the 
theory agree with the 'real' ones describing the simula- 
tion. Generally, a is some function of the control parame- 
ters, a = C[ip]^ and we are looking at the approximate in- 
verse, ipMCTi'p) = CMCT°^i'^~ Ideally, this would 
mean that i^mct = ~ with a = 1. The approxi- 
mate nature of MCT and the approximate matching of 
the particle interactions will induce deviations from this 
ideal behavior. 

Figure W7\ shows the resulting </JMCT(<i5) for the one-, 
three- and five-component analyses we performed. We 
determined this relation by requiring the best possible 
description of the species-averaged collective density cor- 
relators at qd = 7.3 (which essentially involves matching 
the a-rclaxation time scale of the curves). Linear regres- 
sion fits are also shown, and they yield indeed awl, 
so that the fits we shall discuss below are highly rea- 
sonable. Let us point out that we do not see signifi- 
cant deviations from linearity, not even at the ip = 0.585 
that marks the high-density end of our simulations. It 
is usually expected that such deviations set in close to 
ip'^ due to the appearance of hopping processes missed 
in MCT; the commonly reported trends would appear 
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FIG. 13: (Color onhne) MCT fits of the simulated collective 
density correlation functions (shown as circles) using the sim- 
ulated static structure factors binned into M = \ (dashed), 
M = 3 (solid), and M = 5 (dash-dotted lines) compo- 
nents to approximate the simulated polydispersity distribu- 
tion. Packing fractions in the simulation are tp — 0.5, 0.53, 
0.55, 0.57, 0.58, and 0.585. The curves have been fitted by ad- 
justing only v^MCT as described in conjunction with Fig. 1121 
we get ipMCT,M=i = 0.473, 0.502, 0.52, 0.54, 0.554, 0.558; 
V5MCT,Af=3 = 0.449, 0.472, 0.493, 0.5122, 0.5234, 0.5289; and 
V5MCT,M=5 = 0.45, 0.47, 0.49, 0.509, 0.5207, 0.5259. 

in Fig. [TH as a sublinear growth of the (/3MCT-versus- 
(p curve at high densities, owing to the effects of hop- 
ping transport (slower increase in experimental r values 
than predicted by MCT), which would need to be mim- 
icked in MCT by a saturation in the ipuCT variation. 
These trends are typically reported once the a-relaxation 
time exceeds the characteristic time of microscopic short- 
time motion by 10'^ (the most recent claims are for col- 
loidal hard-sphere like systems [13]). Our simulations 
are clearly inside that regime, however we do not find 
such deviations. Note however that many previous stud- 
ies based their conclusions on relaxation times obtained 
for tagged-particle quantities, while we center the dis- 
cussion on collective density correlators. Differences may 
thus partly be attributed to the fact that not all coherent 
incoherent relaxation times diverge in the same manner 
as predicted by MCT. Indeed, fitting the mean-squared 
displacement alone, a different ipmCTif) relation is ob- 
tained [1^ , leading to a different estimate for the critical 
point, and different deviations from linearity. In Fig. 1121 
crosses mark this relation for the M = 3 system. More 
data points are required to address this issue further . 

From the curves shown in Fig. [T^l one can read off the 
estimated critical packing fractions ip"^ for the simulation 
data, as determined from the full MCT analysis once 
Vmct known (the latter are shown as horizontal dotted 
lines). We obtain ip" « 0.594 (0.595, 0.597) for the M = 1 
(A/ = 3; il/ = 5) analysis, i.e. all values agree within less 
than 1% and are in good agreement with the asymptotic 
analysis of the simulation data alone [2^, as expected. 

Figure [Ol shows the simulated coherent correlators for 
qd = 7.3 = qpd together with the MCT curves for the 
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FIG. 14: (Color online) MCT description of the collective den- 
sity correlation functions from the simulation at qd = 12.8 
(corresponding to the second peak of the averaged static 
structure factor). Packing fraction and symbols as in Fig. 1131 



one-component system, as well as the M = 3 and M = 5 
systems, for different densities, from which the relation 
discussed above in conjunction with Fig. [TH has been 
fixed. 

Although little difference can be seen at this wave- 
vector magnitude between the different M, the multi- 
component models give slightly better agreement with 
the simulation, mainly because they show a more 
stretched final decay. As discussed above, the MCT re- 
sults yield Pq values which are too high, but the trend 
with increasing M slightly changes j3q in the right di- 
rection. It is nevertheless remarkable, that the shape of 
the a relaxation is much better reproduced in the MCT 
calculation than one would expect from the difference 
in Pq visible in Fig. [5] (still about 0.2 even at g = qp). 
This clearly indicates that the Kohlrausch function is at 
best an approximate characterization of the a-relaxation 
function at these wave numbers. 

Let us remark that the simulation correlators exhibit 
an unexpected behavior in the very last phase of their de- 
cay at high densities e. g. for Lp = 0.58 and 0.585 at val- 
ues for $(q, t) below 0.05. Here they show a strongly de- 
creasing slope (a 'foot'). The effects on the Kohlrausch- 
fits has been determined to be ±3%, by applying the 
same fit routine and omitting the last part of the curves. 
Thus this behavior cannot be an explanation for the mis- 
match in the stretching exponents. 

Having now fixed all parameters, we compare in Fig. 1141 
simulation and MCT correlators for qd ^ 12.8, corre- 
sponding to the second peak in S{q). The overall fit qual- 
ity is found to be the same as for q — qp, with changes in 
stretching, relaxation time and plateau height that are 
well captured in the theory. Again, curves for different 
M agree, with deviations becoming visible at the high- 
est density investigated. It can also be noted that MCT 
overestimates the correlators in the crossover-region from 
the microscopic to the structural relaxation part. Such 
effects will become more apparent; see below. 
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FIG. 15: (Color online) MCT description of collective den- 
sity correlation functions from the simulation at qd = 9.5 
(corresponding to the first dip after the main peak of the av- 
eraged static structure factor) . Packing fraction and symbols 
as in Fig. 1131 Arrows on the right mark the plateau values 
obtained from Kohlrausch fits discussed in Fig. [5] The two 
curves starting at the plateau values are the solutions to the 
corresponding a-master equation. 



The next point of interest is the first dip of the struc- 
ture factor, qd ~ 9.5. In the region where S{q) < 1 
holds, the differences of the one- and multi-component 
MCT calculations are strongest concerning the plateau 
values (see Section fill C[) . As a consequence of this, the 
time-dependent MCT correlators show stronger discrep- 
ancies for one- and multi-component results. Figure [15] 
shows the comparison of these curves with the simulation 
data for qd ~ 9.5. Indeed, the M = 1 result clearly de- 
viates from the M ~ 3 and M ~ 5 ones. However, none 
of them gives a satisfying description of the simulation 
data in the intermediate time window 0.1 < t < r^; only 
at the longest times, the MCT curves describe the data 
well. This reflects the fact that the g-dependence of is 
well reproduced, see Fig. [71 so that fitting r for g = qp 
also gives good agreement for the final relaxation time at 
other q. 

There is an interesting feature observed in Fig. 1151 
while judging from Fig. [5|the plateau values are in good 
agreement between simulation and MCT for M > 3, this 
agreement is not obvious in the correlators. The reason 
is that in the MCT curves, a pronounced stretched decay 
of the correlators from about (p ~ 0.8 down to (/) « / is 
visible. The simulation curves in contrast show a clear 
shoulder at ss /, following a relatively steep decay from 
the short-time regime. The latter is also well known from 
MD simulations [50| . It is one of the main problems MCT 
has in describing the early /3-relaxation window. 

In general the behavior seen in Figures [TSl and [HI is ex- 
emplary: MCT solutions for the g- values belonging to re- 
gions where S{q) < 1 give worse agreement than the ones 
belonging to g-values where S{q) > 1 holds. The root of 
this problem might be buried in the short-time relax- 
ation part which lasts for longer times at these g-values 
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FIG. 16: (Color online) MOT description of collective density 
correlation functions from the simulation at qd — 3.0, with 
symbols and packing fractions as in Fig. 1131 The inset shows 
only the highest packing fraction, ip — 0.585, and its corre- 
sponding M = 1 correlator, in order to highlight the different 
relaxation times of the different MCT approaches. 



and thus could still influence the /3-process. In one- 
component systems, the short-time relaxation is given 
by ^{q,t) oc exp[—q^/Sqt] for a colloidal system. Hence 
it is conceivable that for lower 5"^, the short-time relax- 
ation is slowed down and hence plays a more important 
role in the dynamics. Since MCT in general only rather 
crudely includes the short-time relaxation, this might be 
the cause of the worse agreement. 

As expected from Fig.Othc multi-component approach 
matches the simulation correlators much better for low 
q values. This is exemplified in Fig. 1161 where data for 
qd = 3.0 are shown. The one-component MCT solution 
underestimates the structural relaxation times by one to 
two decades (cf. the inset of Fig. [TBI), while already M = 3 
gives an overall fit to the data that is much better. This is 
of course another manifestation of the qualitative change 
in Tq for small q discussed in conjunction with Fig. [71 
and hence a signature of the interdiffusion process that 
is absent in the M = 1 calculation. 

Especially at the highest densities, one notices in the 
MCT solutions for M = 3 and M = 5 the emergence of 
a double a-relaxation phenomenon, visible as a shoulder 
around ^ « 0.2. Remarkably, both the M = 3 and i\/ = 5 
results are in close agreement. The simulation data docs 
not show such a double-relaxation pattern, which wc at- 
tribute to the fact that the simulation is truly polydis- 
perse (containing as many species as there are particles), 
and that thus the different a relaxations stemming from 
the superposition of a structural relaxation with different 
interdiffusion processes are smeared out. 

At lower densities, the signature of the interdiffusion 
process remains in the MCT curves as a kink in the relax- 
ation curve for « 0.1, followed by a 'foot' that extends 
almost over one decade in time at the lowest (p shown 
in Fig. 1161 Interestingly, the simulation data for this 
density also show such a feature, and are in fact, apart 
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FIG. 17: (Color onhne) MCT description of the tagged- 
particle density correlation functions determined from the 
simulation at qd = 7.3; packing fractions and symbols as in 
Fig.dSl 

from a shift in time scale for the short-time motion, well 
described by the MCT curve for M > 3. It will be worth- 
while checking for the generality of such a 'foot' in the 
relaxation functions of glass forming systems (as most of 
them are mixtures). 



B. Tagged-Particle Dynamics 

We now address the quality of the MCT description 
for the tagged-particle correlation functions, after all ad- 
justable parameters have been fixed through an analysis 
of the collective density fluctuations at qp. Figure [T71 
shows the results for <f)^{q,t) a.t q ~ qp = 7.3, i.e., it is to 
be compared to Fig.[T21 Two trends are visible in Fig.fTTl 
that mark the main shortcomings of MCT in describ- 
ing the tagged-particle dynamics: first, the relaxation to 
the plateau values is too slow, like in Fig. [151 Secondly, 
a shift in the a-relaxation time is noted, in agreement 
with the expectation from Fig. [71 the MCT curves de- 
cay too slowly, although the collective relaxation times 
match those of the simulation at the same wave number. 
Another general finding for the tagged-particle dynam- 
ics is that the number of component bins M to model 
the polydisperse distribution has little infiuence on the 
quality of the description. 

The agreement in improves somewhat with higher 
q; in agreement with this, also the tagged-particle corre- 
lators are somewhat better described by MCT for larger 
q, as shown for the exemplary case qd — 9.5 in Fig. 1181 
However, the mismatch in the plateau region remains 
roughly the same as in Fig. [T71 

The worsening of the quality of the MCT fits with de- 
creasing q for the tagged-particle quantities was already 
noted in Ref. (29j , and seems to point to an inherent error 
in the theory, that is, however, too poorly understood to 
be improved upon. The error furthermore increases with 
increasing packing fraction, and cannot be eliminated by 
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FIG. 18: (Color online) MCT description of the tagged- 
particle density correlation functions as in Fig. 1171 but for 
qd = 9.5. 




FIG. 19: (Color online) Mean squared displacement from the 
simulation (open circles) compared with their MCT descrip- 
tion based on the fit performed in Fig. 1131 different line styles 
correspond to M = 1, M — 3, and M — 5 models of the 
polydispersity distribution, cf. Fig. 1131 Arrows indicate the 
times for which the van Hove functions are shown in Fig. [20l 
(dashed) and Fig. 1211 (solid arrows). 



alluding to polydispersity effects. It is, of course, di- 
rectly connected to the well known decoupling between 
diffusivity and finite-g local relaxation times, upon which 
we will embark again below. It is therefore not surpris- 
ing that the most drastic deviations between simulation 
and MCT arc visible in the mean-squared displacements, 
shown in Fig. 1191 since this corresponds to the g — ^ 
limit of tagged-particle density correlations. The long- 
time relaxation leading to the final diffusive part in the 
MSD is much slower in MCT than it is in the simulations; 
while at ip = 0.5, for the long time diffusion both curves 
agree within 24%, at = 0.585, the discrepancy is about 
a factor 3. The agreement is again not remedied by in- 
cluding a better description of polydispersity effects; in 
fact the agreement is somewhat worsened in the M = 3 
and M = 5 calculations at (p = 0.585 as compared to the 




FIG. 20: (Color online) P(logSr) = ln(10)47r5r-^G,(r-, t) for 
the times marked by the dashed arrows in Fig. 1191 Open 
circles are the simulation results for ip — 0.5, while lines in- 
dicate the corresponding MCT result for the M = 3 multi- 
component model at ipMCT ~ 0.449: solid black lines denote 
the full distribution, while red-dashed (cyan-short-dashed, 
blue-dashed-dotted) mark the distributions involving only 
small (medium, large) particles. 



one-component analysis. Only the shape of the MSD is 
well described by MCT, again with the caveat that the 
relevant length scale, in this case the height of the plateau 
(indicating the squared cage size up to a trivial prefac- 
tor), is in error in MCT; the theory underestimates the 
localization length by about 10%. This quantitative error 
agrees with the one found for the tagged-particle critical 
amplitude h'^ in Fig. [11] (and cannot be accounted for 
by recalling that the MCT-calculatcd critical density is 
too small) . Fitting only the MSD by adjusting v^mct {f) 
based on this data alone, the simulation data can be de- 
scribed quantitatively (see Ref. [2^ and the discussion 
of Fig. I12p . Similar conclusions have been drawn from 
an analysis of dynamic- light-scattering experiments on 
colloidal suspensions [s^, [E^ ■ 

Following an approach suggested in Refs. [13, [H, [sol - 
Isij we investigated the probability distributions of the 
logarithm of single-particle displacements P{\ogiQ{6r),t) 
at a time t. The appearance of different peaks in 
P{\ogiQ{5r),t) is a result of populations of particles with 
different mobilities, and was suggested as origin of the 
failure of MCT to capture the dynamics of the MSD 
[3 [l^. The probability distribution is directly re- 
lated to the van Hove function via P(logio((5r), i) = 
\n{10)4:irSr^Gs{Sr,t), and its shape is independent of 
time for a Gaussian van Hove function Gs{Sr,t) — 
l/{4TTDt)^/'^exp{~Sr^/4Dt) 

Figure [20] shows P(\ogiQ {6r),t) for the lowest pack- 
ing fraction, ip = 0.5, and its MCT fits by the M = 3 
model, {(pMCT = 0.449). Both simulation and MCT show 
van Hove functions that deviate little from the Gaussian 
expected for ordinary diffusion. As already clear from the 
mean-squared displacements, Fig. I19[ the peak position 
in MCT is at lower Sr values than in the simulation. The 
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FIG. 21: (Color online) P{\og6r) as in Fig. [20] for the times 
marked by solid arrows in Fig. \W\ and for packing fraction 
if — 0.585 (fitted by ipMCT,M=3 = 0.5289). The inset shows 
the average MCT distributions with an additional result at 
t = 1.01 10^ where long-time diffusion has already set in. 
The dotted grey line is a fit with a Gaussian distribution, 
(47rL»t)-^/2 exp[-5rV(4Dt)], where D = 8.35 x 10~® is taken 
from Fig. [19] 
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FIG. 22: (Color online) a-relaxation time scale r determined 
from the master curve of Fig. [3] (red circles), and inverse long- 
time self-diffusion coefficients 1/D extracted from the sim- 
ulated mean-squared displacements (green diamonds). The 
red-solid line shows the power law r ~ |e| ' with 7 = 2.445 
determined from the M = 3 MCT calculations. Inset: 
product Dt evaluated from the simulation data (black up- 
triangles) and from MCT (red down-triangles) scaled by a 
factor to make them comparable. 



MCT calculation is also shown separated into the three 
different particle sizes; quite intuitively, the smaller par- 
ticles are predicted to move faster on average, hence the 
peak visible in Fig. [50] shifts to the right when considering 
only the small particles. One can imagine that the mean- 
squared displacement is a quantity that is dominated by 
motion of fast particles, so that an improvement on the 
theoretical result may be to give stronger weight to their 
displacements. The resulting shift is however seen from 
Fig. [20] to be insufficient to quantitatively explain the 
difference to the simulation data. 

At higher packing fractions, the MCT description 
worsens still, and one starts to see in the simulation 
strong deviations from Gaussian behavior. This is shown 
in Fig.[2l]for the case ip = 0.585 (v^mct = 0.5289). A sec- 
ond shoulder in the distribution P{logiQ{5r),t) at inter- 
mediate times can be interpreted as "hopping-like" mo- 
tion for a certain fraction of particles [13] ■ Such emergent 
two-peak structures are also known from colloidal gels 
[60| , and binary Lennard Jones mixtures [l3, [H, |62| . 

The appearance of dynamical heterogeneities as sig- 
nalled by Fig. [5T] is usually connected with the decou- 
pling of diffusive and collective (viscous) time scales, i.e., 
the breakdown of the Stokes-Einstein (SE) relation men- 
tioned in the introduction. MCT predicts the SE re- 
lation to hold close to ip'^, as both the inverse of the 
long-time self-diffusion coefficient, 1/D, and the typical 
a-relaxation time scale Tq should diverge with the same 
asymptotic power law, so that DTq approaches a con- 
stant as (f ^ (p'^ from below. Figure [H] displays a typ- 
ical collective relaxation time, r, i.e., the one extracted 
from the determination of the a-master curves. Fig. [31 
and the inverse diffusion coefficient as functions of ip for 
the simulation results. Over the window accessible in 



our simulations, both quantities do show power-law-like 
divergence, and in particular r can be reasonably well 
fitted with the expected MCT asymptote, r ^ \o'\~'' , as- 
suming that e cx cr and 7 = 2.445 can be taken from the 
M = 3 calculation. The von Schweidler fits in Fig. [3] 
correspond to 7 « 2.63, so that the MCT exponent re- 
lations are (alomost) fulfilled. The diffusivities, however, 
arc better described by a similar power law with expo- 
nent 7MSD = 2.07942. 

The decoupling of viscous and diffusive time scales is 
best exhibited by taking the product Dt, which is shown 
in the inset of Fig. [22] One clearly notes in the simulation 
data a change of about a factor 4, and no tendency to 
approach a constant at the highest tp we can investigate. 

Rectification plots as shown in Fig.[53]corroborate that 
the asymtotic powerlaw holds for the collective relaxation 
time T over at least two decades, and that a differing 
exponent or a different critical density would be required 
to render the diffusivity as a compatible MCT power law. 

It is instructive to compare the wave-vector depen- 
dence of the relaxation times r" {q) for the tagged-particle 
correlation function with the limiting behavior expected 
on hydrodynamic grounds, T'^{q — >■ 0) ~ 1/q^D. To this 
end, wc plot in Fig.^^ q^T^ (q) extracted from both the 
simulation and our MCT fits, at representative low and 
high packing fractions. Such plots have been suggested 
in discussing the non-Fickian transport evidenced by the 
van Hove functions shown above [l^lsll- There, T'^{q) 
has been defined as the point where the g-dependent 
tagged-particle correlation function has decayed to 1/e. 
This quantity (open symbols and dashed lines in Fig. l24p 
approaches 1 /D at low q. At large q it is expected to drop 
sharply: as the amplitude of the structural relaxation 
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FIG. 23: (Color online) Upper panel; Inverse self Diffusion 
coefficient l/D (green diamonds) and Q-relaxation time scale 
T determined from the master curve of Fig. |3] (red circles) 
with double logarithmic axes. Lower panel: l/D^^''^ (green 
diamonds) and t~^I~< (red circles) from the upper panel. The 
errobars are estimated from the different 7-values obtained 
from the MCT calculations for M = 1,3, 5 components. In 
both panels the solid red line shows the corresponding power 
law from Fig. 1221 




FIG. 24: (Color online) Relaxation times r| of the tagged- 
particle density correlation functions at packing fraction = 
0.585 (upper panel; </3mct = 0.5289) and <^ = 0.5 (lower 
panel; (ySMCT = 0.449), plotted as g^r,. Filled symbols (full 
lines): r| extracted from Kohlrausch fits to the simulation 
(theory) data, taking into account the plateau values /|. 
Open symbols (dotted lines): determined via (^"(g, r^) = 1/e. 



process, drops below 1/e, the procedure no longer re- 
liably probes slow relaxation, but rather is dominated by 
the microscopic short-time relaxation (essentially l/I3o 
in a Brownian system, where l/i'o ^ 1/-D). A simi- 
lar remark holds for the collective relaxation time r((7); 
the time-scale determined by the 1/e-criterion can only 
be compared to MCT as long as /g is sufficiently larger 
than 1/e. 

Our findings for r" (5) are in agreement with the Brow- 



nian dynamics simulations of Flenner and Szamel jlSl j: 
while MCT simply predicts a monotonic crossover be- 
tween the two regimes, in the simulation data, an inter- 
mediate maximum at qd corresponding to the nearest- 
neighbor distance emerges as one approaches . 

Since our MCT fits are matched to the collective cor- 
relation functions, a notable consequence in Fig. [24] is 
that for g — )• 0, MCT and simulation data for q^T^{q) ap- 
proach different constants, the MCT one being too high. 
This corroborates our interpretation that for the tagged- 
particle dynamics MCT is a reasonable theory for the in- 
termediate and large wave numbers, and that deviations 
are increasingly seen as q approaches zero. It differs from 
the interpretation of Ref. |18| , where also simulation and 
MCT data for q^T^{q) were compared, but normahzed 
to their respective g — >■ values. As a result, devia- 
tions were attributed mostly to the intermediate-g, not 
the small-g regime. In light of our results, it might be 
more suggestive to turn around the discussion: it is not 
the nearest-neighbor-scale modes that arc unexpectedly 
slow, but it is the diffusion that is faster than expected 
from the MCT-embedded cage picture. 

If one investigates the distinct part of the collective dy- 
namics, i.e., density correlations that arise from distinct 
particles, MCT's mis-description of tagged-particle dy- 
namics has an interesting consequence. Recall that the 
distinct van Hove correlation function Gd(r,t) is given 
by the difference of collective and tagged-particle contri- 
butions, Gd{r) — G{r) — Ggir)] this means that within 
MCT, it is obtained from the inverse Fourier transform 
of $((/, t) — 4'^{q, t). On physical grounds, Gd{r) must be 
a positive real function, since it measures the probabil- 
ity of finding a particle at time t and distance r from a 
distinct particle that was at the origin at t = 0. This 
property is not obvious from the difference formula. 

Symbols in Fig. [25] show Gd{r,t) evaluated at vari- 
ous times covering the structural-relaxation regime for 
the simulations at the lowest and highest packing frac- 
tion wc studied; we have normalized Gd{r) such that 
it approaches unity at long distances. We recover the 
expected shell structure that is inherited from the n- 
th neighbor shells in the radial distribution function 
g{r) = Gd{r,t = 0). These shells are increasingly 
washed out as time progresses, until the long-time limit 
Gd{r,t ~ 00) = 1 is reached. Comparing with the MCT- 
calculated quantities (obtained from the fits in g-space 
presented above), we note that at distances r including 
and exceeding the nearest-neighbor distance, the MCT 
description is fairly good although not perfect. At small 
r, however, there is a most obvious error as Gd{r, t) turns 
negative in the MCT approximation. Note that for small 
and for sufficiently large i, this phenomenon does not oc- 
cur, for trivial reasons; for small t, the positiveness of 
Gd{r,t) within MCT hinges on that of g{r), while for 
large Gd{r,t) approaches the uniform density. Note 
that the positivity of g{r) may fail for approximate S{q) 
at some densities, but we have checked that this is not 
the case here. 
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FIG. 25: Distinct van Hove correlation fundtions Gdij') at 
packing fractions = 0.5 (upper panel) and 0.585 (lower 
panel) as obtained by simulation (symbols), at various times 
as indicated by arrows in Fig. 1191 Lines are the correspond- 
ing MOT fits evaluated by inverse Fourier transform of the 
difference between collective and tagged-particle density cor- 
relation functions. 



The reason for the failure in the small-r description 
is easily understood: both G{r,t) and Gs{r,t) are dom- 
inated by a strong peak centered on r = 0. since their 
t = values incorporate a (S-peak that is smeared out 
with time. Evaluating Gd{r,t), we have to subtract 
these two large contributions from each other. In fact, 
for the plots shown, typical values of G{r = 0,t) and 
Gs{r = 0,t) are « 14 at the intermediate time, while 
\Gd{r = 0,t)\/G{r = 0,i) = 0(0.1). Numerically, this 
is a moderate error, which by looking at the distinct 
van Hove function is turned into a qualitative one. The 
error was in fact to be expected based on our discus- 
sion so far: the theory, by way of underestimating the 
single-particle diffusion coefficient, overestimates the lo- 
calization of a tagged particle. This translates into a peak 
in Gs {r, t) that is to narrow and thus too high. Even if 
the description of the collective dynamics through G(r, t) 
were totally correct, this overestimation of single-particle 
localization is sufficient to render Gd{r,t) unphysically 




r/d 

FIG. 26: Distinct part of the van Hove correlation func- 
tion Gd{r) that corresponds to the MCT plateau values of 
a monodisperse hard-sphere system within the Percus-Yevick 
approximation, as a function of distance in units of the sphere 
diameters. 



negative. 

In the simulation data, one notices a subtle feature 
around r = that is, due to the reasons just outlined, 
outside the scope of MCT. At r < d, Gd{r, t) raises from 
zero at short times to unity at long times. At the low 
density shown in Fig. I25[ this filling in of the excludcd- 
volume gap happens as naively expected from the broad- 
ening of the nearest-neighbor peak, resulting in functions 
Gd{r,t) that are always monotonically increasing with r 
in the regime r < d. At higher density, however, this 
monotonicity is lost, and an additional dip in Gd{r,t) 
evolves around r = d/2. This qualitatively agrees with 
earlier findin gs f rom simulations of glass-forming binary 
mixtures (6^ |65| . It is intuitively interpreted as the per- 
sistence of preferred interparticle distances: as a given 
particle moves away from its original position, there is 
an enhanced probability that another particle fills this 
position, rather than any nearby one. 

It is worth noting that the issue of negative Gd{r) does 
not appear to be related to the dynamics. As shown 
in Fig. [211 the real-space representation of the noner- 
godicity parameters, calculated from the inverse Fourier 
transform of S{q)f{q) — fiq), shows the same features 
as discussed above. To ensure that no cutoff problems 
arise, we have based this quantity on the Percus-Yevick 
static structure factor with a wave-vector grid spacing 
Aq = 0.05 and M = 1600 discretization points, and 
at packing fraction (p = 0.516 for one-component MCT. 
Note that we are still sufficiently far from the packing 
fraction if > 0.6 where the PY-g(r) starts to be unphys- 
ical. 
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V. CONCLUSION 

We performed molecular-dynamics computer simula- 
tions of a polydisperse quasi-hard sphere system and 
analysed both the collective and the incoherent density- 
density correlation functions in the framework of both 
asymptotic predictions and full numerical solutions of the 
mode-coupling theory of the glass transition. For the lat- 
ter, the required input in the form of static equilibrium 
structure factors has been also calculated from the MD 
simulation. To capture some essential effects of particle- 
size polydispersity in the MCT calculation, numerical so- 
lutions of one-, three-, and five-component systems with 
particle sizes chosen to match the first few moments of 
the true polydispersity distribution have been compared. 

For the particular size distribution chosen in the sim- 
ulation, the five-component analysis turned out to be 
mostly sufficient to capture the effects induced by a vari- 
ation of particle size. These effects concern in particu- 
lar the small-wave-number limit, g — > 0, where mixtures 
of particles of unequal interactions (here: unequal sizes) 
show, even in the species-averaged correlation functions, 
signatures of interdiffusion processes that are well under- 
stood in principle in the framework of multicomponent 
hydrodynamics. MCT has this hydrodynamic limit built 
in, and consequently predicts a subtle interplay of the 
various interdiffusion modes with structural relaxation, 
leading at low q to the appearance of doublc-a peaks. 
These are less pronounced in the simulation, presumably 
due to them being smeared out in a truly polydisperse 
system, where the binning into a few number of compo- 
nents is not strictly meaningful. 

Apart from the low-g behavior, which however may 
be of significant interest as it determines the most com- 
monly discussed transport coefficients, there is little dif- 
ference between the different mixtures we considered. A 
general trend is that, taking the species average already 
on the level of static structure correlations, before enter- 
ing MCT, is worse than performing this average only on 
the dynamic level. This is intuitively understood from 
the fact that such pre-averaging tends to smear out the 
correlations visible in the static structure factor as os- 
cillations at large q. For systems where the slow dy- 
namics is driven by excluded- volume interactions, such 
as our quasi-hard-sphere system, this is a relatively mi- 
nor quantitative effect. For systems where short-range 
interactions are crucial for the slow dynamics (affecting 
the large-q tail of the structure functions), it will be even 
more crucial to treat polydisperse systems as mixtures 
rather than as cffectivc-one-component systems. 

Based on the simulation data and the known asymp- 
totic results from MCT alone, we have first performed 
a traditional scaling analysis, revealing the coefficients 
involved in desribing the slow a relaxation, and the /3- 
relaxation window at intermediate times. The coeffi- 
cients, such as the plateau heights and critical amplitudes 
show good quantitative agreement with the correspond- 
ing quantities calculated within full MCT. A stretched- 



exponential relaxation analysis for the a process reveals 
some mismatch in the stretching indices. 

In fitting the full numerical MCT solutions to the data, 
only the relation between (^mct, the packing fraction 
entering the MCT vertex, and ip, the nominal packing 
fraction in the simulation, was adjusted, to absorb the 
well-known error in the numerical value of ip"^ when cal- 
culated within full MCT. Indeed, the resulting relation 
can be very well described as linear with slope unity, so 
that effectively no fit parameter remains (or just one, if 
one accounts for the slope being slightly different from 
unity) . 

Basing this adjustment on the collective correlation 
functions for a single wave vector magnitude, qd = 7.3, 
we obtain good agreement between MCT and the simu- 
lation data for the collective density correlators at essen- 
tially all sufficiently large qd. Deviations are most promi- 
nent where S{q) < 1. At small qd, the one-component 
analysis shows severe deviations, in particular an order- 
of-magnitude mismatch in the relaxation time, attributed 
to the missing interdiffusion process discussed above. 
The three- and five-component MCT calculations do not 
suffer from this shortcoming and describe also the collec- 
tive small-g behavior reasonably well. 

Wc have discussed the MCT description of the tagged- 
particle dynamics after all parameters have been fixed in 
the analysis of collective correlation functions. Here, the 
situation is more subtle: for qd exceeding roughly half 
the position of the structure-factor main peak, the MCT 
description is again reasonable. However, for g — > 0, 
errors increase continuously, to become most prominent 
in an analysis of the mean-squared displacement (as has 
been noted in an earlier publication). Here, the MCT 
curves show a much stronger slowing down with increas- 
ing packing fraction, while the simulations exhibit av- 
eraged particle mobilities that are higher than those 
expected from either MCT or a Stokes-Einstein argu- 
ment. This diffusion-relaxation decoupling is of course 
well known in the glass literature. The single-particle 
motion has been thoroughly investigated before in terms 
of van Hove functions, identifying subsets of fast and slow 
particles even in one-component systems [6^ . Also in our 
system, such a splitting is observed. However, we wish 
to stress that it seems to affect mostly the MCT descrip- 
tion of tagged-particle dynamics. At the same time, the 
MCT description of the collective density fiuctuations re- 
mains remarkably accurate, and in particular the struc- 
tural relaxation time extracted from the simulation docs 
not show any significant deviation from the values pre- 
dicted by the theory. For future work an analysis of the 
collective van Hove functions might provide additional 
information about the origin of these deviations. 

We can of course not exclude the possibility that at 
even higher densities than those we could simulate, such 
deviations eventually set in. However, if this is the case, 
they need not coincide with the features typically dis- 
cussed in terms of heterogeneous dynamics, viz. the de- 
coupling of diffusivity from structural relaxation. It may 
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be that our choice of a system driven by stochastic dy- 
namics and without a significant energy scale in the parti- 
cle interactions is fortuitous. However, this remains to be 
clarified. Within MCT, the independence of structural- 
relaxation properties of the time-evolution operator is a 
major result, confirmed before [sO]. There are arguments 
that this correspondence does indeed extend also to the 
way simulation results deviate from MCT, but not to 
higher-order correlation functions such as four-point sus- 
ceptibilities [63, HI. 

It appears then that fitting incoherent correlation func- 
tions with MCT is a rather roundabout way of testing the 
theory, and that in particular it represents an unfortu- 
nate test (for the theory, at least) in that these quantities 
show strongest deviations from the predicted behavior. 
Unfortunately, the MSD is perhaps the quantity most 
often analyzed in simulation and in assessing MCT's va- 
lidity (since it is easy to compute with good statistics). 
Our analysis shows that it is the least valuable quantity 
to compare MCT with. 

Our discussion suggests that other collective correla- 



tion functions provide a fruitful basis for further inves- 
tigations of MCT and its approximations: in particu- 
lar stress-stress auto-correlation functions should be an- 
alyzed, as they are cloesly linked to the memory kernel 
of the theory. However, these require even more com- 
putational effort to determine from simulation than the 
collective density correlators we have examined here: for 
the latter, we have averaged 1000 independent evalua- 
tions, while for the former, up to 4000 had to be used 
in previous work on a similar system with Newtonian 
dynamics [69j . 
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